Расчет статически неопределимой стержневой системы

Условие

>    restart;

>    n:=[1,2,4]:

>    with(LinearAlgebra):

Координаты шарнира балки

>    x0:=0: y0:=0:

Координаты шарниров на балке

>    xb:=<40,20,60>: yb:=<10,5,15>:

Координаты опорных шарниров

>    xe:=<0,20,20>: ye:=<10,45,-25>:

Длины стержней

>    for i to 3 do

>     v[i]:=Vector([xe[i]-xb[i],ye[i]-yb[i],0]):

>     L[i]:=VectorNorm(v[i],2);

>     v1[i]:=s[i]*Normalize(v[i],Euclidean);

>     r[i]:=Vector([xb[i]-x0,yb[i]-y0,0]):

>    od:

Координаты точки приложения силы

>    r[4]:=Vector([80,20,0]):

Сила

>    v1[4]:=Vector([16,0,0]):

 Метод сил

>    N:=Matrix(3):

>    for j to 3 do

>      i:=n[j]:s[i]:=1:

>      EQ:=CrossProduct(r[3],v1[3])+CrossProduct(r[i],v1[i]):

>      N[j,j]:=1: N[j,3]:=evalf(solve(EQ[3],s[3]));

>    od:

Интеграл Мора

>    for i to 2 do

>    for j to 3 do

>    d[i,j]:=evalf(add(L[k]*N[i,k]*N[j,k],k=1..3));od;od;

Каноническая система метода сил

>    for i to 2 do

>    eq[i]:=d[i,1]*x1+d[i,2]*x2+d[i,3];od:

>    sl:=solve({eq[1],eq[2]},{x1,x2}):assign(sl);

>    s[1]:=x1:s[2]:=x2:

Уравнение моментов

>    EQ:=add(CrossProduct(r[k],v1[k]),k=1..4):

>    x3:=evalf(solve(EQ[3],s[3])):